!>\brief
!! Darcy test case
!!
!! We solve a simple diffusion equation where the diffusion
!! coefficient has a high variability in space, in order to simulate
!! the presence of permeable and impermeable materials. The two
!! diffusion coefficients are \f$\varepsilon_1\f$ and
!! \f$\varepsilon_2\f$, with 
!! \f{displaymath}{
!!  \frac{\varepsilon_2}{\varepsilon_1}\ll 1.
!! \f}
!<----------------------------------------------------------------------
module mod_darcy_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_darcy_test_constructor, &
   mod_darcy_test_destructor,  &
   mod_darcy_test_initialized, &
   test_name, &
   test_description,&
   coeff_diff,&
   coeff_adv, &
   coeff_re,  &
   coeff_f,   &
   coeff_dir, &
   coeff_neu, &
   coeff_rob

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "darcy"

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(2)
 logical, protected ::               &
   mod_darcy_test_initialized = .false.
 ! private members
 real(wp), parameter :: &
   eps1 = 1.0_wp, &
   eps2 = 1e-6_wp
 integer :: ndim
 character(len=*), parameter :: &
   this_mod_name = 'mod_darcy_test'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_darcy_test_constructor(d)
  integer, intent(in) :: d

  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_darcy_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ndim = d
   write(test_description(1),'(a)') &
     'Darcy test case'
   write(test_description(2),'(a,e9.3,a,e9.3,a)') &
     '  viscosities eps1 = ',eps1,', eps2 = ',eps2,'.'

   mod_darcy_test_initialized = .true.
 end subroutine mod_darcy_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_darcy_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_darcy_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_darcy_test_initialized = .false.
 end subroutine mod_darcy_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_diff(x) result(mu)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: mu(size(x,1),size(x,1),size(x,2))

  real(wp), parameter :: u = 0.1_wp*(5.0_wp/4.0_wp), & ! unit
    u1=1.0_wp*u, u2=2.0_wp*u, u3=3.0_wp*u, &
    u5=5.0_wp*u, u6=6.0_wp*u, u7=7.0_wp*u

  integer :: l, i
  real(wp) :: xx, yy, zz, eps
 
   mu = 0.0_wp
   do l=1,size(x,2)

     xx = x(1,l)
     yy = x(2,l)

     select case(ndim)

      case(2)
       if( ((yy>0.6_wp).and.(yy<0.7_wp).and.(xx<0.9_wp)) .or. &
           ((yy>0.3_wp).and.(yy<0.4_wp).and.(xx>0.1_wp)) ) then
         eps = eps2
       else
         eps = eps1
       endif

      case(3)
       zz = x(3,l)
       if( ((zz>u1).and.(zz<u3).and.((xx>u2).or.(yy>u1))) .or. &
           ((zz>u5).and.(zz<u7).and.((xx<u7).or.(yy<u6))) ) then
         eps = eps2
       else
         eps = eps1
       endif

     end select

     do i=1,size(x,1)
       mu(i,i,l) = eps
     enddo

   enddo
       
 end function coeff_diff
 
!-----------------------------------------------------------------------
 
 pure function coeff_adv(x) result(a)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: a(size(x,1),size(x,2))

   a = 0.0_wp
 end function coeff_adv
 
!-----------------------------------------------------------------------

 pure function coeff_re(x) result(sigma)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: sigma(size(x,2))
 
   sigma = 0.0_wp
 end function coeff_re
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,2))

   f = 0.0_wp
 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,2))

   bregcase: select case(breg)
    case(1)
     d = 1.0_wp
    case(3,6) ! 2D and 3D
     d = 0.0_wp
    case default
     d = huge(x(1,1))
   end select bregcase
 
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_neu(x,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: h(size(x,2))

   bregcase: select case(breg)
    case(2,3,4,5)
     h = 0.0_wp
    case default
     h = huge(x(1,1))
   end select bregcase
 
 end function coeff_neu
 
!-----------------------------------------------------------------------

 pure function coeff_rob(x,breg) result(g)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: g(size(x,2))

   g = 0.0_wp
 end function coeff_rob
 
!-----------------------------------------------------------------------

end module mod_darcy_test

